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Abstract 

Large-scale  simulations  are  often  under-resolved  at  some  level,  but 
they  are  still  useful  in  extracting  both  qualitative  and  quantitative  infor¬ 
mation  about  the  flow.  In  order  to  use  such  results  effectively  we  need 
to  characterize  the  numerical  uncertainty  of  under-resolved  simulations. 
However,  different  numerical  methods  exhibit  different  behavior,  and 
spectral-based  methods  in  particular  may  over-predict  fluctuations  both 
in  amplitude  and  frequency  due  to  their  very  low  artificial  dissipation 
in  contrast  with  finite  differences.  In  this  chapter,  we  provide  insight 
into  under-resolved  spectral  simulations  and  document  several  diagnos¬ 
tic  signs  of  under-resolution  for  spectral/hp  element  methods.  We  first 
review  the  state-of-the  art  in  direct  numerical  simulation  and  present 
a  new  class  of  spectral  methods  on  unstructured  grids  for  handling 
complex-geometry  compressible  and  incompressible  flows.  We  then  fo¬ 
cus  on  the  effects  of  under-resolving  the  nonlinear  contributions,  and 
finally  we  present  prototype  cases  for  both  transitional  and  turbulent 
flows. 


Keywords:  Spectral  methods,  complex-geometry,  under-resolution,  unstructured 
grids,  turbulence 
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1.  Introduction 

Under-resolved  simulations  are  perhaps  the  rule  rather  than  the  excep¬ 
tion!  This  can  be  understood,  as  in  practice  users  attempt  high  Reynolds 
number  simulations  in  problems  with  new  physics  and  thus  unknown 
resolution  requirements.  Verification  and  validation  of  the  solution  is  a 
very  tedious  process  ([1]),  and  at  present  there  are  no  established  efficient 
methods  to  assess  numerical  accuracy.  Also,  for  large-scale  turbulence 
simulations,  existing  computational  resources  may  often  be  inadequate 
for  the  attempted  simulations,  so  additional  error  checking  simulations 
would  be  prohibitively  expensive. 

However,  an  under-resolved  simulation  is  not  useless  but,  in  fact,  it 
can  provide  a  lot  of  information  if  proper  characterization  is  established 
combined  with  experience  for  the  specific  discretization  used  in  such  a 
simulation.  One  such  example  is  a  relatively  early  direct  numerical  sim¬ 
ulation  of  turbulent  channel  flow  by  ([2])  which  has  remained  largely 
unnoticed.  In  figure  1  we  plot  one  of  his  results,  i.e.  the  Reynolds  stress 
distribution  across  the  channel.  It  is  in  good  agreement  with  high  resolu¬ 
tion  DNS  even  for  the  lowest  resolution  employed  in  Zores’s  simulation, 
corresponding  to  four  Fourier  modes  in  the  streamwise,  16  Chebyshev 
modes  in  the  normal,  and  six  Fourier  modes  in  the  spanwise  direction. 
The  Reynolds  number  based  on  the  wall  shear  velocity  is  R*  ~  120.  In 
order  to  achieve  smooth  profiles  a  very  long-  time  averaging  was  em¬ 
ployed.  Clearly,  this  is  an  example  of  an  under-resolved  simulation, 
which  however  sustains  the  turbulence  fluctuations  and  leads  to  better 
than  10%  accuracy  in  second-order  statistics.  In  contrast,  a  low  reso¬ 
lution  simulation  based  on  finite  differences  would  typically  converge  to 
the  laminar  steady  state  solution. 

For  under-resolved  simulations  to  be  useful  we  need  to  characterize 
both  numerical  and  physical  uncertainty,  creating  appropriate  composite 
error  bars  similar  to  experiments.  This  is  a  very  difficult  task,  and  work 
on  uncertainty  associated  with  the  data,  i.e.  input,  is  still  at  an  early 
stage.  On  the  numerical  side,  there  are  still  many  classical  issues  which 
are  unresolved  today,  e.g.  skew-symmetry  of  advection  operators  in  the 
discrete  sense,  time-integration  algorithms  with  large  time  step,  efficient 
treatment  of  geometric  complexity,  efficient  adaptivity,  etc. 

There  are  two  major  challenges  today  in  direct  numerical  simula¬ 
tions  (DNS)  of  turbulence  following  the  successes  of  the  last  two  decades 
( [3] ) :  The  first  is  that  the  maximum  Reynolds  number  possible  in  simula¬ 
tions  is  still  much  lower  compared  to  turbulent  flows  of  practical  interest. 
For  example,  at  present  or  in  the  near  future,  the  maximum  Re \  (based 
on  the  Taylor  micro-scale)  for  homogeneous  turbulence  that  can  be  ac- 
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curately  simulated  is  less  than  500.  However,  in  geophysical  flows  the 
typical  Reynolds  number  Re\  may  be  orders  of  magnitude  higher.  The 
second  challenge  we  face  is  that  complex-geometry  flows  are  still  largely 
untackled;  geometries  beyond  the  standard  channel  flow  with  flat  walls 
have  only  recently  been  considered,  and  most  of  them  involve  at  least 
one  homogeneous  direction. 

A  summary  of  the  range  of  Reynolds  number  and  geometries  for  which 
direct  numerical  simulations  have  been  successfully  completed  is  plotted 
in  the  sketch  of  figure  2.  It  shows  that  accurate  direct  numerical  simula¬ 
tions  of  turbulence  in  simple-geometry  domains  can  handle  much  higher 
Reynolds  number  flows  than  in  complex-geometry  domains.  This,  in 
essence,  reflects  the  additional  computational  complexity  associated  with 
discretization  of  the  Navier-Stokes  equations  in  complex-geometry  do¬ 
mains.  Clearly,  the  Fourier  discretization  which  is  employed  for  all  three 
directions  in  homogeneous  turbulence  cannot  be  used  in  inhomogeneous 
directions,  where  Chebyshev  or  Legendre  spectral  discretizations  (or  al¬ 
ternatively  some  high-order  finite  difference  variant)  are  used.  More 
specifically,  on  non-separable  or  multiply-connected  domains  (e.g.  flow 
past  a  circular  cylinder)  these  classical  methods  are  also  inappropriate, 
and  thus  new  domain-decomposition  based  methods  need  to  be  used 
effectively. 

As  regards  the  type  of  discretization,  it  was  evident  even  from  the 
early  attempts  to  perform  DNS  of  turbulence  in  the  seventies,  that 
high-order  discretization  was  not  only  computationally  advantageous  but 


Figure  1.  Reynolds  stress  distribution  of  a  low-resolution  simulation  of  a  turbulent 
channel  flow  by  ([2]).  The  different  symbols  correspond  to  different  resolution  in  x 
(stream),  y  (normal)  and  z  (span)  as  follows:  A  :  8 x 32  x 8;  ixj:  8 x  16 x 8;  +  :  6 x  16 x  6; 
x  :  4  x  16  x  6. 
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Geometric  Complexity 

Figure  2.  Conceptual  overview  of  DNS  of  turbulent  flows:  Maximum  Reynolds 
number  versus  geometric  complexity. 

also  a  necessity.  Simulating  turbulence  requires  long-time  integration, 
however  non-negligible  dispersion  errors  associated  with  low-order  dis¬ 
cretization  could  eventually  render  the  computational  results  erroneous. 
There  is  plenty  of  anecdotal  evidence  about  such  results  from  the  early 
practitioners,  and  modern  numerical  analysis  can  rigorously  document 
this  as  well.  The  importance  of  high-order  discretization  has  also  been 
recognized  for  large  eddy  simulations  (LES)  ([4]),  as  discretization  errors 
seem  to  interact  with  the  subgrid  modeling  errors  in  an  adverse  way. 

As  we  go  through  the  multi- Teraflop  (and  beyond)  computing  era 
and  are  capable  of  performing  simulations  of  1  billion  points  or  more 
at  reasonable  turn-around  time,  high-order  numerical  methods  will  play 
a  key  role  in  simulating  high  Reynolds  number  and  complex-geometry 
turbulence.  They  provide 

■  fast  convergence, 

■  small  diffusion  and  dispersion  errors, 

■  easier  implementation  of  the  inf-sup  condition  for  incompressible 
Navier-Stokes, 
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■  better  data  volume-over-surface  ratio  for  efficient  parallel  process¬ 
ing,  and 

■  better  input /output  handling  due  to  the  smaller  volume  of  data. 

For  many  engineering  applications  where  accuracy  of  the  order  of  10% 
is  acceptable,  quadratic  convergence  is  usually  sufficient  for  stationary 
problems.  However,  this  is  not  true  in  time-dependent  flow  simulations 
where  long-time  integration  is  required.  Also,  in  DNS  a  10%  inaccuracy 
in  phase  errors  may  lead  to  flow  re-laminarization.  Therefore,  we  must 
ask  how  long-time  integration  relates  to  the  formal  order  of  accuracy  of 
a  numerical  scheme,  and  what  is  the  corresponding  computational  cost? 
To  this  end,  let  us  consider  the  convection  of  a  waveform  at  a  constant 
speed.  Let  us  now  assume  that  there  are  grid  points  required  per 
wavelength  to  reduce  the  error  to  a  level  e,  where  k  denotes  the  formal 
order  of  the  scheme.  In  addition,  let  us  assume  that  we  integrate  for  M 
time  periods.  We  can  neglect  temporal  errors  0(At)J  (where  J  is  the 
order  of  the  time  integration)  by  assuming  a  sufficiently  small  time  step 
At.  We  wish  to  estimate  the  phase  error  in  this  simulation  for  second- 
N^2\  fourth-  IV(4),  and  sixth-  order  finite  difference  schemes.  The 
complete  analysis  is  presented  in  ([5]),  and  here  we  present  the  results 
for  the  computational  work.  In  figure  3  we  compare  the  efficiency  of 
these  three  different  discretizations  for  the  same  phase  error  by  plotting 
the  computational  work  required  to  maintain  an  “engineering”  accuracy 
of  10%  versus  the  number  of  time  periods  for  the  integration.  This 
comparison  favors  the  fourth-order  scheme  for  short  times  ( M  oc  0(1)) 
over  both  the  second-order  and  the  sixth-order  schemes.  However,  for 
long-time  integration  (M  oc  0(100)),  even  for  this  engineering  accuracy 
of  10%,  the  sixth-order  scheme  is  superior  as  the  corresponding  operation 
count  H%6)  is  about  6  times  lower  than  the  operation  count  of  the  second- 
order  scheme  W^2\  and  half  the  work  of  the  fourth-order  scheme  % 
For  an  accuracy  of  1%  in  the  solution  of  this  convection  problem,  the 
sixth-order  scheme  is  superior  even  for  short-time  integration. 

High-order  accuracy,  however,  does  not  automatically  imply  a  resolved 
and  thus  accurate  DNS  or  LES.  In  particular,  spectral-based  methods 
tend  to  behave  differently  when  the  number  of  grid  points  or  modes  is 
insufficient.  For  example,  they  tend  to  be  more  unstable,  lead  to  over¬ 
prediction  of  amplitudes,  and  could  even  result  in  erroneous  unsteady 
flow  at  subcritical  conditions.  This  is  the  primary  topic  that  we  focus 
on  in  the  present  paper. 

Specifically,  we  first  review  some  key  developments  in  extending  spec¬ 
tral  methods  to  unstructured  grids  for  both  incompressible  and  com¬ 
pressible  flows.  We  then  discuss  in  some  detail  the  effect  of  under- 
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Figure  3.  Computational  work  (number  of  floating-point  operations)  required  to 
integrate  a  linear  advection  equation  for  M  periods  while  maintaining  a  cumulative 
phase  error  of  10%. 


resolving  the  discretization  of  nonlinear  terms  and  how  dealising  can  be 
handled  on  non-uniform  grids  presenting  both  one-dinrensional  exam¬ 
ples  but  also  full  DNS  of  turbulent  flow  that  may  suffer  from  aliasing. 
We  then  proceed  with  several  examples  of  internal  and  external  flows 
and  document  diagnostics  that  can  be  employed  to  detect  erroneous 
physics.  We  also  include  simulations  of  some  turbulent  flows  which,  al¬ 
though  clearly  under-resolved,  lead  to  useful  results  in  agreement  with 
the  experiments.  Finally,  we  conclude  with  a  perspective  on  simulat¬ 
ing  turbulence  in  fully  three-dimensional  domains  where  both  numerical 
uncertainty  and  physical  uncertainty  are  adequately  characterized. 

2.  Spectral  Methods  on  Unstructured  Grids 

There  have  been  more  than  fifteen  years  of  developments  in  extending 
spectral  methods  to  complex-geometry  domains  ([5]),  starting  with  the 
pioneering  work  of  ([6]),  who  developed  spectral  methods  in  the  con¬ 
text  of  a  multi-element  variational  formulation  similar  to  finite  element 
methods.  This  allowed  the  use  of  spectral  (Chebyshev  or  Legendre) 
expansions  as  trial  basis  in  general  quadrilateral  subdomains.  Conti- 
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nuity  of  data  and  unknowns  across  subdomains  is  ensured  via  appro¬ 
priate  construction  of  the  trial  basis  similar  to  finite  element  methods 
for  second-order  elliptic  problems.  Such  methods  were  used  to  produce 
the  first  spectral  DNS  of  turbulence  in  complex-geometry  domain,  flow 
over  riblets,  in  ([7]).  An  extension  to  non-conforming  discretizations  for 
turbulent  flows,  which  are  more  appropriate  for  local  refinement,  was 
presented  in  ([8]). 

The  new  generation  of  spectral  methods  developed  recently  is  more 
appropriate  for  discretizations  on  unstructured  grids  consisting  of  tri¬ 
angles  and  tetrahedra,  similar  to  grids  used  in  aerodynamics  ([9,  10]). 
In  many  simulations,  however,  it  is  more  efficient  to  employ  hybrid  dis¬ 
cretizations,  i.e.  discretizations  using  a  combination  of  structured  and 
unstructured  subdomains.  This  is  a  recent  trend  in  computational  me¬ 
chanics  involving  complex  three-dimensional  computational  domains  ( 
[11,  12]).  Such  an  approach  combines  the  simplicity  and  convenience 
of  structured  domains  with  the  geometric  flexibility  of  an  unstructured 
discretization.  In  two-dimensions,  hybrid  discretization  simply  implies 
the  use  of  triangular  and  rectangular  subdomains,  however  in  three- 
dimensions  the  hybrid  strategy  is  more  complex  requiring  the  use  of 
hexahedra,  prisms,  pyramids  and  tetrahedra. 

We  have  developed  a  unified  description  in  dealing  with  elements  of 
different  shape  in  two-  and  three-dimensions.  This  unified  approach 
generates  polynomial  expansions  which  can  be  expressed  in  terms  of  a 
generalized  product  of  the  form 

<j)pqr(x,y,z)  =  (j)ap{x)4>hpq{y)4>cpqr{z). 

Here  we  have  used  the  Cartesian  co-ordinates  x,  y  and  z  but,  in  general, 
they  can  be  any  set  of  co-ordinates  defining  a  specified  region.  The  stan¬ 
dard  tensor  product  is  simply  a  degenerate  case  of  this  product,  where 
the  second  and  third  functions  are  only  dependent  on  one  index.  The 
primary  motivation  in  developing  an  expansion  of  this  form  is  computa¬ 
tional  efficiency.  Such  expansions  can  be  evaluated  in  three-dimensions 
in  0(P4)  operations  as  compared  to  0(P 6)  operations  with  non-tensor 
products  (where  P  is  the  number  of  spectral  modes  per  direction). 

2.1  Local  Co-ordinate  Systems 

We  start  by  defining  a  convenient  set  of  local  co-ordinates  upon  which 
we  can  construct  the  expansions.  Unlike  the  barycentric  co-ordinates, 
which  are  typically  applied  to  unstructured  domains  in  linear  finite 
elements,  we  define  a  set  of  collapsed  Cartesian  co-ordinates  in  non- 
rectangular  domains.  These  co-ordinates  will  form  the  foundation  of  the 
polynomial  expansions.  The  advantage  of  this  system  is  that  every  do- 


Figure  4-  Triangle  to  rectangle  transformation 


Figure  5.  Hexahedron  to  tetrahedron  transformation 


main  can  be  bounded  by  constant  limits  of  the  new  local  co-ordinates; 
accordingly  operations  such  as  integration  and  differentiation  can  be 
performed  using  standard  one-dinrensional  techniques. 

The  new  co-ordinate  systems  are  based  upon  the  transformation  of 
a  triangular  region  to  a  rectangular  domain  (and  vice  versa)  as  shown 
in  figure  4.  The  main  effect  of  the  transformation  is  to  map  the  verti¬ 
cal  lines  in  the  rectangular  domain  (i.e.  lines  of  constant  rj\)  onto  lines 
radiating  out  of  the  point  (£i  =  —1,^2  =  1)  in  the  triangular  domain. 
The  triangular  region  can  now  be  described  using  the  “ray”  co-ordinate 
(?/i)  and  the  standard  horizontal  co-ordinate  (£2  =  772).  The  triangu¬ 
lar  domain  is  therefore  defined  by  (—1  <  771,7/2  <  1)  rather  than  the 
Cartesian  description  (—1  <  £1,  £2;  £1  +  ^2  <  0)  where  the  upper  bound 
couples  the  two  co-ordinates.  The  “ray”  co-ordinate  (77 1)  is  multi-valued 
at  (£1  =  — 1,^2  =  !)•  Nevertheless,  we  note  that  the  use  of  singu¬ 
lar  co-ordinate  systems  is  very  common,  arising  in  both  cylindrical  and 
spherical  co-ordinate  systems. 

As  illustrated  in  figure  5,  the  same  transformation  can  be  repeatedly 
applied  to  generate  new  co-ordinate  systems  in  three-dimensions.  Here, 
we  start  from  the  bi-unit  hexahedral  domain  and  apply  the  triangle  to 
rectangle  transformation  in  the  vertical  plane  to  generate  a  prismatic 
region.  The  transformation  is  then  used  in  the  second  vertical  plane  to 
generate  the  pyramidic  region.  Finally,  the  rectangle  to  triangle  trans- 
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Figure  6.  The  local  coordinate  systems  used  in  each  of  the  hybrid  elements  and  the 
planes  described  by  fixing  each  local  co-ordinate. 

formation  is  applied  to  every  square  cross  section  parallel  to  the  base  of 
the  pyramidic  region  to  arrive  at  the  tetrahedral  domain. 

By  determining  the  hexahedral  co-ordinates  (r/i .  7/2,  r/3)  in  terms  of 
the  Cartesian  co-ordinates  of  the  tetrahedral  region  (£1 ,  £2,  £3)  we  can 
generate  a  new  co-ordinate  system  for  the  tetrahedron.  This  new  system 
and  the  planes  described  by  fixing  the  local  co-ordinates  are  shown  in 
figure  6.  Also  shown  are  the  new  systems  for  the  intermediate  domains 
which  are  generated  in  the  same  fashion.  Here  we  have  assumed  that 
the  local  Cartesian  co-ordinates  for  every  domain  are  (£i,£2>£3)- 

2.2  Hierarchical  Expansions 

For  each  of  the  hybrid  domains  we  can  develop  a  polynomial  expansion 
based  upon  the  local  co-ordinate  system  derived  in  section  2.1.  These 
expansions  will  be  polynomials  in  terms  of  the  local  co-ordinates  as  well 
as  the  Cartesian  co-ordinates  (^1,^2,  £3)-  This  is  a  significant  property 
as  primary  operations  such  as  integration  and  differentiation  can  be 
performed  with  respect  to  the  local  co-ordinates  but  the  expansion  may 
still  be  considered  as  a  polynomial  expansion  in  terms  of  the  Cartesian 
system. 

We  shall  initially  consider  expansions  which  are  orthogonal  in  the 
Legendre  inner  product.  We  define  three  principle  functions  ^“(2:),  fiiji2) 
and  4>ijp(z),  in  terms  of  the  Jacobi  polynomial,  P“’/3(z),  as: 

=  r°'%),  4>%F)  =  (t)‘  p? w'V). 

=  (A) ’+3  g2,+2J+2,0(T 
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Using  these  functions  we  can  construct  the  orthogonal  polynomial  ex¬ 
pansions: 

Hexahedral  expansion:  <?W(£  1,62,  £3)  =  (6)  <^(6  )</£(£  3) 

Prismatic  expansion:  0pgr(6>6>6 3)  =  </>p(6)<(,q(r7  2)<^r(6) 
Pyramidic  expansion:  0pgr(6>6>6)  =  <t>p(W)&aq(m)<t>CpqArlz) 
Tetrahedral  expansion:  0P9r(6,  6> 6)  = 
where, 


»7i 


2(1+6) 

(-6  -  6) 


m 


2(1  +  6) 
(1  -  6) 


% 


2(1  +  6) 
(1  -  6) 


%  =  6, 


are  the  local  co-ordinates  illustrated  in  figure  6. 

The  hexahedral  expansion  is  simply  a  standard  tensor  product  of 
Legendre  polynomials  (since  Pp,0(z)  =  Lp(z)).  In  the  other  expansions 
the  introduction  of  the  degenerate  local  co-ordinate  systems  is  linked  to 
the  use  of  the  more  unusual  functions  4>\j(z)  and  <j>°jk(z).  These  func¬ 
tions  both  contain  factors  of  the  form  j  which  is  necessary  to  keep 
the  expansion  as  a  polynomial  of  the  Cartesian  co-ordinates  (6, 6, 6)- 
For  example,  the  co-ordinate  r) 2  in  the  prismatic  expansion  necessitates 
the  use  of  the  function  </>j)r(6)  which  introduces  a  factor  of  ^ 1 j  .  The 
product  of  this  factor  with  ^>“(7/2)  is  a  polynomial  function  in  6  and  6- 
Since  the  remaining  part  of  the  prismatic  expansion,  </>“( 6),  is  already 
in  terms  of  a  Cartesian  co-ordinate  the  whole  expansion  is  a  polynomial 
in  terms  of  the  Cartesian  system. 

The  polynomial  space,  in  Cartesian  co-ordinates,  for  each  expansion 
is: 

V  =  Span{£?  $  6}  (1) 

where  pqr  for  each  domain  is 


Hexahedron 

Prism 

Pyramidic 

Tetrahedron 


0  <  p  <  Pi 
0  <  p  <  Pi 
0  <  p  <  Pi 
0  <  p  <  Pi 


0<q<P2 
0  <q<P2 
0<q<P2 
0  <  p  +  q  <  +2 


0  <  r  <  P3 

0  <  q  +  r  <  P3  . 

0<p  +  q  +  r<Ps  ^ 

0  <  p  +  q  +  r  <  P3. 


The  range  of  the  p,  q  and  r  indices  indicate  how  the  expansions  should 
be  expanded  to  generate  a  complete  polynomial  space.  We  note  that  if 
P\  =  P2  =  +3  then  the  tetrahedral  and  pyramidic  expansions  span  the 
same  space  and  are  in  a  subspace  of  the  prismatic  expansion  which  is  in 
turn  a  subspace  of  the  hexahedral  expansion. 
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2.3  Galerkin  and  Discontinuous  Galerkin 
Projections 

To  obtain  a  system  of  nonlinear  algebraic  equations,  we  employ  dif¬ 
ferent  projections  and  time  integration  algorithms.  In  particular,  for 
incompressible  flows  we  use  a  linear  Galerkin  projection  in  conjunction 
with  the  high-order  fractional  stepping  scheme  described  in  ([13,  14]). 
For  compressible  flows,  we  use  a  discontinuous  Galerkin  projection  with 
multi-step  explicit  time  integration  ([15]). 

We  describe  both  approaches  next  with  more  emphasis  on  the  latter 
which  is  a  more  recent  development. 


Incompressible  Flows.  The  standard  approach  in  treating  the 
incompressible  Navier-Stokes  equations  is  to  combine  a  semi-implicit 
scheme  with  a  fractional  procedure  ([5])  following  the  Eulerian  descrip¬ 
tion.  Here  we  consider  a  more  recent  development  that  takes  advantage 
of  semi-Lagrangian  treatment  for  advection.  This  allows  for  large  size 
time  steps  in  simulations  of  turbulence,  where  at  high  Reynolds  num¬ 
ber  the  temporal  scales  are  largely  over-resolved.  Following  ([16])  we 
consider  the  Navier-Stokes  equations  in  Lagrangian  form 

du  _  _ •)  .  . 

—  =  -Vp  +  i/V“  u,  (3) 


V  •  u  =  0,  (4) 

where  d/dtt  denotes  a  Lagrangian  derivative.  We  employ  a  stiffly-stable 
second-order  scheme  to  discretize  the  time  derivative: 


|  un+1  -  2 u'i  +  M"1 


=  (-Vp  +  zzV2«y 


where  u 
and 


i 


is  the  velocity  u  at  the  departure  point  x%  at  time  level  tn, 
is  the  velocity  at  the  departure  point  x^~l  at  time  level  fn_1. 


The  departure  point  x %  is  obtained  by  solving 


—  =un+l/2{x,t),  x(tn+1)  =  xt 

at 

and  also 

un+l/2  =  3/2un  _  1/2 un~l 

The  point  a^-1  is  obtained  by  solving 

^  =  un(x,t),  x(tn+1 )  =  xa. 
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By  using  the  above  characteristic  equations,  the  resulting  scheme  is 
second-order  accurate  in  time. 

Specifically,  for  computational  convenience  we  use  the  following  three 


substeps  to  solve  equation  (5) 

fi-2«3  +  i«r1  n 

At 

(6) 

*rfi  = 

At. 

(7) 

—  —  11 

U  =  uV2un+1. 

At 

(8) 

The  discrete  divergence-free  condition  results  in  a  consistent  Poisson 
equation  for  the  pressure,  i.e. 


V  V+1  =  -J-V  •  u, 

1  At 

with  accurate  pressure  boundary  conditions  of  the  form  ([5]) 

on 

where  n  is  the  unit  normal,  and  oj  is  the  vorticity  ([14]). 

The  semi-Lagrangian  approach  is  typically  more  expensive  than  the 
corresponding  Eulerian  approach,  but  in  practice  the  larger  size  of  time 
step  allowed  in  the  former  leads  to  more  efficient  simulations.  This  was 
demonstrated  for  two-  and  three-dimensional  flows  in  ([16]). 

With  regards  to  spatial  discretization ,  in  order  to  enforce  the  required 
C°  continuity,  the  orthogonal  expansion  is  modified  by  decomposing 
the  expansion  into  an  interior  and  boundary  contribution.  This  results 
in  partially  sacrificing  orthogonality.  The  interior  modes  (or  bubble 
functions)  are  defined  to  be  zero  on  the  boundary  of  the  local  domain. 
The  completeness  of  the  expansion  is  then  ensured  by  adding  boundary 
modes  which  consist  of 

■  Vertex,  Edge,  and  Face  contributions. 

The  vertex  modes  have  unit  value  at  one  vertex  and  decay  to  zero  at 
all  other  vertices;  edge  modes  have  local  support  along  one  edge  and  are 
zero  on  all  other  edges  and  vertices,  and  face  modes  have  local  support 
on  one  face  and  are  zero  on  all  other  faces,  edges  and  vertices.  C° 
continuity  between  elements  can  then  be  enforced  by  matching  similar 
shaped  boundary  modes.  The  local  co-ordinate  systems  do  impose  some 
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restrictions  on  the  orientation  in  which  triangular  faces  may  connect. 
However,  it  has  been  shown  in  ([10])  that  a  C°  tetrahedral  expansion 
can  be  constructed  for  any  tetrahedral  mesh.  A  similar  strategy  could 
be  applied  to  a  hybrid  discretization. 

Compressible  Flows.  We  consider  the  non-dimensionalized  com¬ 
pressible  Navier-Stokes  equations,  which  we  write  in  compact  form  as 

Ut  +  V  •  F  =  R(QlV  ■  ¥u  (9) 

where  F  and  Fu  correspond  to  inviscid  and  viscous  flux  contributions, 
respectively.  Here  the  vector  U  =  [p,  pu,  pv,  piu,  E]1  with  ( u,v,w )  the 
local  fluid  velocity,  p  the  fluid  density,  and  E  the  total  internal  en¬ 
ergy.  Splitting  the  Navier-Stokes  operator  in  this  form  allows  for  the 
separate  treatment  of  the  inviscid  and  viscous  contributions,  which  in 
general  exhibit  different  mathematical  properties.  In  the  following,  we 
review  briefly  the  discontinuous  Galerkin  formulations  employed  in  the 
proposed  method.  A  systematic  analysis  of  the  advection  operator  was 
presented  in  ([17]),  where  a  mixed  formulation  was  used  to  treat  the 
diffusion  terms.  No  flux  limiters  are  necessary  as  has  been  found  before 
in  ([18])  and  has  been  justified  theoretically  in  ([19]). 

We  first  use  a  linear  two-dimensional  advection  equation  of  a  con¬ 
served  quantity  u  in  a  region  H,  in  order  to  illustrate  the  treatment  of 
inviscid  contributions: 


—  +  V-F(«)  =  0,  (10) 

where  F (u)  =  ( f(u),g(u),h(u ))  is  the  flux  vector  which  defines  the 
transport  of  u(x,f).  We  start  with  the  variational  statement  of  the 
standard  Galerkin  formulation  of  (10)  by  multiplying  by  a  test  function 
v  and  integrating  by  parts 

[  ^vdx+  [  vh-F(u)ds  —  f  Vv-F(u)dx  =  0.  (11) 

Jq  ot  Jen  Jn 

The  solution  u  G  X  (approximation  space)  satisfies  this  equation  for 
all  v  £  V  (test  space).  The  requirement  that  X  consist  of  continuous 
functions  naturally  leads  to  a  basis  consisting  of  functions  with  over¬ 
lapping  support,  which  implies  equation  (11)  becomes  a  banded  matrix 
equation.  Solving  the  corresponding  large  system  is  not  a  trivial  task  for 
parallel  implementations,  and  therefore  a  different  type  of  formulation 
is  desirable. 
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Figure  7.  Interface  conditions  between  two  adjacent  triangles. 

Another  consideration  from  the  point  of  view  of  advection  is  that  con¬ 
tinuous  function  spaces  are  not  the  natural  place  to  pose  the  problem. 
Mathematically,  hyperbolic  problems  of  this  type  tend  to  have  solutions 
in  spaces  of  bounded  variation.  In  physical  problems,  the  best  one  can 
hope  for  in  practice  is  that  solutions  will  be  piecewise  continuous,  that 
is,  be  smooth  in  regions  separated  by  discontinuities  (shocks).  An  ad¬ 
ditional  consideration  is  that  the  formulation  presented  next  preserves 
automatically  conservativity  in  the  element-wise  sense. 

These  considerations  suggest  immediately  a  formulation  where  X  may 
contain  discontinuous  functions.  The  discrete  space  Xs  contains  poly¬ 
nomials  within  each  “element,”  but  zero  outside  the  element.  Here  the 
“element”  is,  for  example,  an  individual  triangular  region  T%  in  the  com¬ 
putational  mesh  applied  to  the  problem.  Thus  the  computational  do¬ 
main  fl  =  Uj  Tj,  and  T%,  Tj  overlap  only  on  edges. 

Contending  with  the  discontinuities  requires  a  somewhat  different  ap¬ 
proach  to  the  variational  formulation.  Each  element  ( E )  is  treated  sep¬ 
arately,  giving  a  variational  statement  (after  integrating  by  parts  once 
more) : 

St(u,v)e  +  [  v(f(ui,ue)  -  f(ui))-nds  +  (V  •  f(u),v)E  =  0,  (12) 

at  Jote 

where  f{ui)  is  the  flux  of  the  interior  values.  Computations  on  each 
element  are  performed  separately,  and  the  connection  between  elements 
is  a  result  of  the  way  boundary  conditions  are  applied.  Here,  bound¬ 
ary  conditions  are  enforced  via  the  numerical  surface  flux  f(v,i,ue)  that 
appears  in  equation  (12).  Because  this  value  is  computed  at  the  bound¬ 
ary  between  adjacent  elements,  it  may  be  computed  from  the  value  of  u 
given  at  either  element.  These  two  possible  values  are  denoted  here  as 
Ui  in  the  interior  of  the  element  under  consideration  and  and  ue  in  the 
exterior  (see  figure  7).  Upwinding  considerations  dictate  how  this  flux 
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is  computed.  In  the  more  complicated  case  of  a  hyperbolic  system  of 
equations,  an  approximate  Riemann  solver  should  be  used  to  compute 
a  value  of  f,g,h  (in  three-dimensions)  based  on  u,{  and  ue.  Specifically, 
we  compute  the  flux  f(ui,ue)  using  upwinding,  i.e. 

f(u)  =  RA+Lui  +  RA~  Lue 

where  A  (the  Jacobian  matrix  of  F )  is  written  in  terms  of  the  left  and 
right  eigenvectors,  i.e.  A  =  RAL  with  A  containing  the  corresponding 
eigenvalues  in  the  diagonal;  also,  =  (A  ±  |A|)/2.  Alternatively,  we 
can  use  a  standard  Lax-Friedrichs  flux 

f(u)  =  ^(/K)  +  -  ^R\A\L(ue  -  Ui). 

This  last  form  is  what  is  used  in  the  airfoil  example  presented  in  section 
4. 


Next,  we  consider  as  a  model  problem  the  parabolic  equation  with 
variable  coefficient  v  to  demonstrate  the  treatment  of  the  viscous  con¬ 
tributions: 

ut  =  V  •  {u\7u)  +  /,  in  11,  u  £  L2(Q) 
u  =  g(x,  t),  on  dUl 
We  then  introduce  the  flux  variable 

q  =  —vVu 

with  q(x,  t)  £  L2(fl),  and  re-write  the  parabolic  equation 

Ut  =  —V  •  q  +  /,  in  O 

1/i/q  =  —  Vu,  in  Ul 
u  =  g(x,  t),  on  dUl. 

The  weak  formulation  of  the  problem  is  then  as  follows:  Find  (q,  u)  £ 
L2(fi)  x  L2(Ul)  such  that 

(ut,  w)E  =  (q,  Vrc)E-  <  w,  qb  •  n  >E  +(/,  w)E,  Vrc  £  L2(Ul) 

l/z/(qm,  v)E  =  (u,  V  •  v)E-  <  ub,  v  •  n  >E,  Vv  £  L2(f2) 
u  =  p(x,  t),  on  dUl 

where  the  parentheses  denote  standard  inner  product  in  an  element  (E) 
and  the  angle  brackets  denote  boundary  terms  on  each  element,  with  n 
denoting  the  unit  outwards  normal.  The  surface  terms  contain  weighted 
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boundary  values  of  Vb,qb,  which  can  be  chosen  as  the  arithmetic  mean 
of  values  from  the  two  sides  of  the  boundary,  i.e. 

Vb  =  (vi  +  ve)/2, 


and 

qb  =  ( qi  +  <?e)/2. 

The  consequences  of  choosing  different  numerical  fluxes  with  regards  to 
stability  and  accuracy  have  been  investigated  in  ([20]). 

By  integrating  by  parts  once  more,  we  obtain  an  equivalent  formula¬ 
tion  which  is  easier  to  implement,  and  it  is  actually  used  in  the  computer 
code.  The  new  variational  problem  is 

(ut,  w)E  =  (— V  •  q,  w)E~  <  w,  (qb  -  qi)  •  n  >E  +(/,  w)E,  Vrc  G  L2(fl) 

l/v>(q,v)E  =  (-Vu,v)£-  <  ub  -  Ui,  v  •  n  >E,  Vv  G  L2(fi) 
u  =  g(x,  t ),  in  dQ 

where  the  subscript  (i)  denotes  contributions  evaluated  at  the  interior 
side  of  the  boundary.  We  integrate  the  above  system  explicitly  in  time 
and  employ  the  orthogonal  Jacobi  polynomials  as  trial  and  test  basis. 

3.  Nonlinearities  and  Dealising 

In  spectral  methods  the  quadratic  nonlinearities  in  the  incompressible 
Navier-Stokes  equations  or  the  cubic  nonlinearities  in  the  compressible 
Navier-Stokes  are  computed  in  the  physical  space.  Specifically,  the  fields 
(velocity,  pressure,  energy)  are  first  transformed  into  physical  space  and 
subsequently  the  products  are  obtained  at  all  quadrature  points  in  a 
collocation  fashion.  Another  transform  is  then  performed  to  bring  the 
results  back  to  modal  space.  More  specifically,  when  the  number  of 
quadrature  points  Q  is  the  same  as  the  number  of  modes  in  the  spectral 
expansion  P  we  have  a  true  collocation  method,  otherwise  for  Q  >  P 
we  have  a  super-collocation  method. 

The  form  in  which  we  write  the  nonlinear  terms,  that  is, 

■  in  convective  (flux),  or 

■  skew-symmetric,  or 

■  rotation  form 

is  also  important.  In  spectral  DNS  of  boundary  layers  and  channel  flows, 
the  rotation  form  is  usually  preferred  over  the  convective  form  as  it  semi¬ 
conserves  energy  (in  the  inviscid  limit).  This,  typically,  makes  it  more 
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stable,  especially  for  the  long-time  integration  required  in  DNS.  In  ad¬ 
dition,  it  is  more  economical  as  it  requires  the  evaluation  of  only  six 
derivatives  whereas  the  convective  form  requires  nine  derivative  evalu¬ 
ations.  The  skew-symmetric  form  was  found  to  be  more  “forgiving”  in 
aliasing  errors  in  under-resolved  simulations  of  homogeneous  turbulence 
compared  to  the  rotation  form.  This  is  also  true  for  finite  difference 
methods;  see  the  article  of  ([21])  in  this  volume  where  it  is  shown  that 
skew-symmetry  leads  to  symmetry  preservation  and  enhanced  stability. 
However,  the  skew-symmetric  form  requires  the  evaluation  of  18  deriva¬ 
tives  which  is  computationally  more  expensive. 

There  is  not  sufficient  experience  yet  with  spectral/ hp  element  DNS 
to  conclusively  suggest  one  form  or  the  other  although  there  is  some 
consensus  that  the  convective  form  is  quite  accurate  and  leads  to  sta¬ 
ble  discretizations.  A  comparison  of  the  different  forms  (convective, 
fiux,  and  skew-symmetric)  was  performed  in  ([22])  for  a  constant  ad- 
vection  velocity  as  well  as  for  a  spatially  varying  divergent-free  veloc¬ 
ity  ( u,v )  =  (— sin.X2  cos  .Ti, sin xi  cos  .T2).  The  discretization  was  based 
on  a  nodal  Gauss-Lobatto-Legendre  basis.  The  result  was  that  for  the 
constant  advection  velocity  all  forms  were  the  same  in  that  they  pro¬ 
duced  identical  eigenspectrum  with  all  imaginary  eigenvalues.  However, 
for  the  variable  advection  velocity,  only  the  skew-symmetric  form  gave 
imaginary  eigenvalues  with  the  convective  and  conservative  form  produc¬ 
ing  complex  eigenvalues  with  positive  real  parts.  For  purely  convection 
equations,  these  spurious  positive  values  may  lead  to  instabilities  if  ex¬ 
plicit  time  stepping  is  used.  However,  for  Navier-Stokes  computations  at 
modest  Reynolds  number  no  such  instabilities  have  been  observed,  pre¬ 
sumably  due  to  the  stabilizing  role  of  the  viscous  terms.  Although  the 
skew-symmetric  form  is  usually  considered  the  most  accurate,  problems 
may  also  be  encountered  with  this  form  for  Dirichlet  and  inflow/outflow 
conditions.  This  is  presented  in  some  detail  in  ([5])  using  a  numerical 
experiment  performed  by  ([23]). 

Finally,  errors  may  be  caused  by  insufficient  quadrature  used  in  the 
spectral//ip  element  discretization  of  the  nonlinear  terms,  especially  in 
complex-geometry  flows.  These  errors  can  be  eliminated  effectively  by 
employing  over-integration,  i.e.  integrating  the  nonlinear  terms  in  the 
variational  statement  with  higher  order  quadrature  than  the  one  em¬ 
ployed  for  the  linear  contributions,  e.g.  pressure  and  viscous  terms.  We 
will  examine  this  issue  in  some  detail  next. 
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3.1  Accuracy,  Stability  and  Over-Integration 

To  understand  the  ramifications  of  under-integration  of  nonlinear  terms, 
we  perform  the  following  test: 

1  Consider  a  single  element  in  the  space  interval  [—1,1]  containing 
P  =16  Jacobi  modes. 

2  Initialize  all  the  modal  coefficients  to  one. 

3  Evaluate  the  modal  representation  on  a  set  of  Q  quadrature  points. 

4  Square  (in  a  pointwise  fashion)  the  values  at  the  quadrature  points. 

5  Pre-multiply  the  set  of  points  (as  a  vector)  by  the  collocation  deriv¬ 
ative  matrix  of  the  appropriate  size  (rank  Q  x  Q). 

6  Project  back  to  modal  coefficients  by  discrete  inner  products  using 
Gaussian  integration. 

The  procedure  above  mimics  the  “physical  space”  or  pseudo-spectral 
evaluation  of  the  term  commonly  used  in  spectral  methods  for  evalu¬ 
ating  nonlinear  terms.  This  test  was  chosen  because  even  in  its  simplicity 
it  models  the  order  of  nonlinearity  that  occurs  in  the  solution  of  the  in¬ 
compressible  Navier-Stokes  equations.  All  modes  are  set  to  one  to  mimic 
a  case  in  which  an  element  has  under-resolved  or  marginally  resolved  the 
solution  within  the  element.  In  the  test  above,  the  only  unspecified  pa¬ 
rameter  is  the  number  of  quadrature  points  Q  to  be  used.  In  using 
Gauss-Lobattto  points,  the  value  of  Q  is  taken  to  be  one  more  than  the 
number  of  modes  P  (in  this  case  then  P  =  16  and  Q  =  17)  ([24]),  but 
this  value  is  appropriate  for  the  inner  products  corresponding  to  linear 
terms.  For  quadratic  or  cubic  nonlinearities  more  quadrature  points  are 
required.  The  ramifications  of  under-integration  of  this  form  are  shown 
in  figure  8.  The  figure  on  the  left  was  obtained  for  quadratic  nonlinearity 
(■ 2)  and  the  figure  on  the  right  was  obtained  for  a  cubic  nonlinearity 
(■ TfeU 3).  The  difference  in  the  modal  coefficients  at  the  conclusion  of  the 
algorithm  above  for  different  values  of  Q  is  provided.  We  observe  that 
for  the  quadratic  nonlinearity,  once  \P  quadrature  points  are  used,  the 
differences  in  the  modal  values  do  not  change.  Similarly  for  the  cubic 
nonlinearity,  once  2 P  quadrature  points  are  used  the  differences  in  the 
modal  values  do  not  change. 

In  order  to  appreciate  the  effect  of  under-integration  in  the  context  of 
a  numerical  solution,  we  consider  the  inviscid  Burgers  equation,  which 
we  discretize  using  the  discontinuous  Galerkin  method.  The  initial  con¬ 
dition  is  —  sin(7rx),  and  five  equally  spaced  elements  spanning  [—1,1] 
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x  104 


Figure  8.  Comparison  of  the  difference  in  modal  coefficients  when  different  numbers 
of  quadrature  points  are  used.  Quadratic  nonlinearity  is  shown  on  the  left  and  cubic 
nonlinearity  is  shown  on  the  right. 


were  used,  each  one  having  P  =  16  modes.  In  figure  9,  we  plot  the 
L2  norm  of  the  solution  versus  the  number  of  quadrature  points  used 
for  numerical  integration.  When  using  Q  =  17, 19  and  Q  =  21  points, 
the  solution  is  unstable  (denoted  by  the  blue  *).  Once  the  number  of 
quadrature  points  reaches  Q  =  24  (|P  where  P  is  the  number  of  modes), 
the  L2  norm  of  the  solution  does  not  change. 

We  can  analyze  this  behavior  by  examining  the  energy  in  the  modes 
(denoted  by  the  square  of  the  modal  values)  within  the  element  that 
contains  the  jump  in  the  inviscid  Burgers  solution.  The  modes  were 
extracted  at  time  T  =  0.35,  after  the  shock  has  formed  (at  time  ^)  and 
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Figure  9.  Solution  of  the  inviscid  Burgers  equation  evaluated  at  T  =  0.5.  Five 
equal  spaced  elements  were  used  with  16  modes  in  each  element.  On  the  ordinate 
we  plot  the  L2  norm  of  the  solution,  and  on  the  abscissa  we  plot  the  number  of 
quadrature  points  used  for  numerical  integration.  Unstable  solutions  are  denoted  by 
blue  *.  Observe  that  after  Q  =  24  points,  the  L2  norm  of  the  solution  does  not 
change. 


prior  to  the  solution  becoming  unstable.  In  figure  10,  we  plot  the  square 
of  the  modal  coefficients  versus  the  mode  number.  Due  to  the  symmetry 
of  the  element  placement,  only  even  number  modes  were  excited. 

This  case  corresponds  to  Q  =  17  quadrature  points  being  used,  which 
from  the  figure  9,  we  know  will  become  unstable  by  time  T  =  0.5.  If  a 
| P  rule  is  used,  yielding  Q  =  24  points,  the  solution  is  stable,  and  the 
energy  is  much  less  than  when  the  non-linear  terms  are  under-integrated. 
This  plot  shows  vividly  the  effects  of  aliasing  when  under-integration  of 
the  non-linear  terms  is  performed. 

An  alternative  way  of  handling  instabilities  associated  with  nonlinear¬ 
ities  in  hyperbolic  conservation  laws  is  through  monotonicity  preserving 
schemes.  An  approach  suitable  for  high-order  methods  has  been  devel¬ 
oped  by  ([25])  and  was  employed  in  large-eddy  simulations  in  ([26]).  It 
involves  the  addition  of  a  second-order  convolution  kernel  that  acts  on 
each  mode  separately  and  controls  the  high  modes  suppressing  prefer¬ 
entially  erroneous  high-frequency  oscillations.  This  type  of  nonlinear 
kernel  has  been  termed  as  spectral  vanishing  viscosity  (SVV).  As  an 
example,  the  inviscid  Burgers  equation  with  the  SVV  term  added  is 


du  1  du1  9  du 

~dt+ 


(13) 
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Figure  10.  Modal  coefficients  of  the  inviscid  Burgers  solution  before  blow  up.  Both 
over-integration  and  SVV  lead  to  a  stable  solution  unlike  the  collocation  approach. 


Here  e  oc  1/P,  i.e.  it  is  inversely  proportional  to  the  number  of  modes, 
and  Qk  is  a  smooth  kernel  that  facilitates  a  transition  between  the  con¬ 
trolled  high  modes  and  the  uncontrolled  and  more  energetic  low  modes. 
It  is  given  by 

(fc-p)2 

Qk  =  e  <k-pcU  ,  k  >  Pc. 

The  cut-off  wave  number  Pc  scales  as  Pc  ~  \[P,  asymptotically  for  P 
large.  Although  the  above  can  be  considered  as  a  viscosity  regularization 
procedure  there  is  a  significant  difference  as  has  been  demonstrated  in 
([27]),  see  also  ([26]).  The  parameters  e  and  Qk  are  chosen  so  that 
monotonicity  is  preserved  while  the  spectral  accuracy  in  the  solution  is 
also  maintained.  In  figure  10  we  demonstrate  how  the  addition  of  SVV 
leads  to  effectively  the  same  results  as  over-integration  but  operating 
with  a  collocation  discretization,  i.e.  Q  =  P  +  1.  The  modal  coefficients 
of  the  inviscid  Burgers  solution  converge  monotonically  to  zero  leading 
to  a  stable  simulation  unlike  the  collocation  untreated  simulation. 

3.2  Transition  and  Turbulence  in  a  Triangular 
Duct 

We  demonstrate  next  the  effect  of  under-integration  and  associated 
aliasing  errors  by  simulating  transition  to  turbulence  of  incompressible 
flow  in  a  duct  with  its  cross-section  being  an  equilateral  triangle.  The 
laminar  fully-developed  solution  is  known  analytically.  We  introduce 
some  random  disturbances  in  the  flow  and  we  integrate  in  time  until 
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these  disturbances  start  decaying  or  growing  in  time.  All  simulations 
were  performed  in  the  domain  shown  in  figure  11  with  the  cross-section 
discretized  using  one  triangular  element  only  and  16  Fourier  modes  (32 
collocation  points)  in  the  streamwise  (homogeneous)  direction.  The 
Reynolds  number  is  defined  as  Re  =  UDe/v  where  U  is  the  average 
velocity  and  De  is  the  equivalent  (hydraulic)  diameter.  For  Re  <  500  all 
disturbances  decay  but  for  Re  =  1250  the  flow  goes  through  transition 
and  a  turbulent  state  is  sustained. 

We  have  performed  three  simulations  at  Re  =  1250  corresponding 
to  three  different  combinations  of  polynomial  and  quadrature  order.  In 
the  first  one,  shown  in  11(a),  we  consider  the  case  where  Q  =  P  +  1, 
where  P  =  16.  The  forces  on  the  three  walls  of  the  duct  are  plotted 
as  a  function  of  time.  From  symmetry  considerations,  we  expect  that 
the  statistical  averages  of  the  three  forces  are  identical  but  obviously 
the  symmetry  is  not  preserved  here.  In  figure  11(b)  we  plot  the  forces 
for  the  case  with  Q  =  2 P,  and  in  figure  11(c)  the  case  with  Q  =  3P/2. 
We  have  verified  that  in  both  cases  the  same  statistical  force  average 
is  obtained,  consistent  with  the  analysis  presented  above  for  handling 
under- integration  induced  errors. 

Based  on  the  above  analysis  and  result  as  well  as  other  similar  results, 
we  can  state  the  following  semi-empirical  rule: 

Dealising  Rule:  For  quadratic  nonlinearities  employing  super-collocation 
with  3/2  P  grid  (quadrature)  points  per  direction,  where  P  is  the  poly¬ 
nomial  order  per  direction,  followed  by  a  Galerkin  projection  leads  to  a 
dealiased  turbulence  simulation  on  non-uniform  meshes. 

4.  Under- Resolution  and  Diagnostics 

Spectral  and  spectral//ip  element  methods  behave,  in  general,  differ¬ 
ently  than  low-order  methods  in  under-resolved  simulations.  Spectral 
discretizations  are  more  susceptible  to  numerical  instabilities  than  other 
low-order  discretizations.  This  could  be  frustrating  for  the  users  who 
seek  robustness  but  it  is  actually  safe-guarding  against  erroneous  an¬ 
swers,  as  typically  spectral  codes  blow  up  in  seriously  under-resolved 
simulations.  In  under-resolved  spectral  discretizations  there  are  many 
more  wiggles,  and  therefore  it  is  easier  to  detect  suspicious  simulations 
before  resorting  to  more  rigorous  error  estimation.  Also,  spectral  dis¬ 
cretizations  typically  suffer  from  little  numerical  dissipation  unlike  finite- 
difference  methods,  which  introduce  an  erroneous  numerical  viscosity 
in  low  resolution  discretizations.  This  effectively  lowers  the  nominal 
Reynolds  number  of  the  simulated  flow  and  leads  to  stable  simulations 
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Figure  11.  Duct  flow  domain:  The  cross-section  is  an  equilateral  triangle  and  the 
streamwise  length  is  three  times  the  triangle  edge.  Shown  is  a  snapshot  of  streamwise 
velocity  contours  at  Re  =  1250. 


but  with  the  incorrect  physics.  This  is  not  true  in  spectral/hp  discretiza¬ 
tions  where  the  nominal  Reynolds  number  is  also  the  effective  Reynolds 
number.  However,  such  behavior  in  conjunction  also  with  the  presence 
of  high- wave  number  wiggles,  may  sometimes  be  the  source  of  erroneous 
instabilities  in  under-resolved  spectral  flow  simulations.  For  example, 
the  resulting  velocity  profiles  may  not  be  monotonic  and  thus  are  sus¬ 
ceptible  to  inviscid  type  instabilities,  which  in  turn  promote  transition 
from  steady  to  unsteady  flow  or  transition  to  higher  bifurcations  and 
eventually  turbulence.  For  open  unsteady  flows,  the  amplitude  of  the 
oscillation  in  an  under-resolved  simulation  is  usually  over- predicted. 

In  the  following,  we  present  a  few  examples  of  under-resolved  flows 
that  affect  both  transition  to  turbulence  as  well  as  turbulence  statistics. 
We  also  include  a  case  of  transient  compressible  flow  past  an  airfoil  where 
under-resolution  may  seriously  affect  the  lift.  However,  as  we  will  see  not 
all  results  from  under-resolved  simulations  are  inaccurate.  Some  flows 
exhibit  low-dimensionality  and  the  energetics  of  low  modes  dominate  so 
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Figure  12.  Wall  shear  forces  on  each  wall  as  a  function  of  time  for  (a)  (Q  =  M  +  1) 
(b)  (Q  =  2 M);  and  (c)  (Q  =  3M/2). 
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even  a  very  coarse  grid  turbulence  simulation  may  predict  the  correct 
statistics  -  this  is  the  case  of  the  cylinder  wake. 


Figure  13.  Low  resolution  mesh  for  flow  over  a  backwards-facing  step. 
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Figure  14 ■  High  resolution  mesh  for  flow  over  a  backwards-facing  step. 
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Figure  15.  Time  history  at  the  third  point  (shown  in  the  mesh  in  figures  13  and 
14)  at  Re=700.  Solid  line:  low  resolution;  Dash  line:  high  resolution. 

4.1  Erroneous  Flow  Transition 

The  first  example  is  from  the  systematic  spectral  simulations  pre¬ 
sented  by  ([28]),  in  the  study  of  bypass  transition  in  a  boundary  layer. 
Using  a  Chebyshev  discretization  in  the  inhomogeneous  direction  and 
Fourier  expansions  in  the  other  two  directions,  they  demonstrated  that 
with  P  =  33  Chebyshev  modes  their  simulation  showed  that  a  wave 
structure  was  present  which  would  give  rise  to  a  secondary  instabil¬ 
ity,  as  suggested  earlier  by  other  investigators.  However,  this  structure 
changed,  and  the  instability  vanished  completely  when  P  =  66  Cheby¬ 
shev  modes  were  employed  in  the  simulation.  In  fact,  simulations  with 
even  higher  resolution  confirmed  this  explanation. 

An  example  of  similar  behavior  but  with  spectral  element  discretiza¬ 
tion  is  the  simulation  of  flow  over  a  backwards-facing  step,  which  was 
first  presented  in  ([29]).  For  the  resolution  shown  in  figure  13  the  flow  is 
unsteady  as  is  evident  in  the  plot  that  shows  time  history  of  velocity  (fig¬ 
ure  15).  However,  if  higher  resolution  is  used  as  shown  in  figure  14,  then 
a  steady  state  is  reached  (see  corresponding  figure  15),  and  this  is  true 
even  at  higher  Reynolds  numbers  up  to  about  Re  ~  2, 500.  Interestingly, 
the  results  of  the  under-resolved  simulation  are  not  totally  irrelevant  as 
they  contain  information  about  the  actual  flow  albeit  at  a  different  set 
of  parameters.  For  example,  the  two  frequencies  present  in  the  unsteady 
case  are  the  natural  frequencies  of  the  flow  corresponding  to  the  shear 
layer  instability  at  the  step  corner  and  the  Tollmien-Schlichting  waves 
in  the  downstream  channel  portion.  These  modes  are  excited  either  by 
background  noise,  for  example,  some  small  turbulence  level  at  the  in- 


Under-Resolution  and  Diagnostics  in  Spectral  Simulations 


27 


flow,  or  spontaneously  at  a  higher  Reynolds  number.  Since  no  absolutely 
quiet  wind  tunnels  exist,  the  results  of  the  under-resolved  “noisy”  simu¬ 
lation,  in  this  case,  match  the  results  of  the  experiment  ([30]).  We  note 
that  other  inherently  noisy  discretizations  employing  vortex  methods 
and  lattice-Boltzmann  methods  also  lead  to  an  unsteady  flow  solution  ( 
[31,  32]). 

4.2  Fast-Pitching  Airfoil 

Next,  we  consider  laminar  flow  around  a  rapidly  pitching  airfoil  and 
compare  discontinuous  Galerkin  spectral//ip  element  results  against  the 
finite  volume  results  obtained  in  ([33]).  In  particular,  we  consider  a 
NACA  0015  airfoil  pitching  upwards  about  a  fixed  axis  at  a  constant 
rate  from  zero  incidence  to  a  maximum  angle  of  attack  of  approximately 
60  degrees.  The  pivot  axis  location  is  at  1/4  of  the  chord  measured  from 
the  leading  edge.  The  temporal  variation  of  the  pitch  given  in  ([33])  is 

n(t)  =  n0[ i -e"4-6t/to],  t>o 

where  to  denotes  the  time  elapsed  for  the  airfoil  to  reach  99%  of  its  final 
pitch  rate  fig.  Here  the  non-dimensional  values  are  fg  =  1.0  and  Hg  =  0.6 
based  on  the  chord  length  and  free  stream  velocity.  As  initial  condition 
the  computed  field  at  0  degrees  angle  of  attack  is  used.  The  Mach 
number  is  M  =  0.2  and  the  chord  Reynolds  number  is  Re  =  10, 000. 

In  ([33])  a  similar  simulation  was  obtained  using  a  grid  fixed  to  the 
airfoil  by  employing  an  appropriate  transformation  and  discretizing  the 
modified  compressible  Navier-Stokes  equations  using  the  implicit  ap¬ 
proximate  factorization  of  ([34]).  A  typical  grid  used  in  ([33])  involved 
203  x  101  points.  In  the  present  study,  we  employ  the  domain  shown  in 
figure  16.  We  performed  two  different  sets  of  simulations,  first  with  un¬ 
structured  discretization  around  the  airfoil  (see  figure  17;  total  of  3,888 
triangular  elements),  and  subsequently  with  hybrid  discretization  with 
quadrilateral  elements  around  the  airfoil  for  better  resolution  of  bound¬ 
ary  layers  (total  of  116  quadrilateral  and  2167  triangular  elements).  We 
demonstrate  how  the  hybrid  discretization  combined  with  variable  P- 
order  per  element  allows  accurate  resolution  of  boundary  elements  with¬ 
out  the  need  for  re-meshing.  We  first  performed  simulations  with  con¬ 
stant  P-order  on  all  elements  and  subsequently  with  higher  P-order  in 
the  inner  layers  of  elements  as  shown  in  figure  18.  We  contrast  the  results 
in  figure  19  for  P-order  P  =  3  on  the  left,  and  P  varying  from  10  in  the 
innermost  layer  to  2  in  the  far  field.  We  see  that  the  boundary  layer  is 
unresolved  as  indicated  by  the  discontinuities  at  the  element  interfaces, 
but  it  is  accurately  resolved  in  the  second  simulation. 


Figure  18.  Hybrid  discretization  showing  the  variable  p-order  on  a  gray-scale  map 
around  the  airfoil.  The  dash  vertical  line  indicates  the  location  where  boundary  layer 
profiles  are  taken  (see  figure  19). 
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Figure  19.  Boundary  layer  profiles  for  a  simulation  with  uniform  P-resolution  (left) 
and  variable  P-resolution  (right,  as  shown  in  figure  18). 
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Figure  20.  Lift  (upper  curve)  and  drag  (lower  curve)  coefficients  versus  angle  of 
attack  in  degrees.  The  symbols  correspond  to  computations  of  ([33]),  the  dot  fine 
corresponds  to  our  simulation  at  P  =  2,  the  solid  fine  to  P  =  3  and  the  dash  line  to 
P  =  4. 


Returning  now  to  the  unstructured  grid,  we  test  convergence  by  also 
performing  R-refinenrent  on  the  same  triangulization  but  with  three  dif¬ 
ferent  values  of  spectral  order  R  corresponding  to  2nd,  3rd  and  4th  order 
polynomial  interpolation.  In  figure  20  we  plot  the  computed  lift  and  drag 
coefficients  versus  the  angle  of  attack  for  grids  corresponding  to  P  =  2, 3 
and  P  =  4.  We  also  include  (with  symbols)  the  computational  results 
of  ([33]),  and  we  see  that  in  general  there  is  very  good  agreement  except 
at  the  large  angles  of  attack  close  to  50  degrees.  This  difference  is  due 
to  qualitative  difference  in  ffow  structure  at  small  scales,  which  are  only 
resolved  with  the  higher  order  simulations. 

To  examine  differences  in  the  ffow  field  due  to  spatial  resolution  we 
plot  in  figure  21  density  contours  for  the  cases  P  =  2  and  P  =  3  at 
non-dimensional  time  t  =  0.75  corresponding  to  an  angle  of  attack  18.55 
degrees.  We  see  that  the  higher  resolution  simulation  provides  a  more 
detailed  picture  of  the  vortex  shedding  in  the  near-wake,  but  the  contours 
around  the  airfoil  are  very  similar.  At  a  later  time  t  =  1.5,  corresponding 
to  an  angle  of  attack  of  44.1  degrees,  there  are  differences  between  the 
computations  at  resolution  P  =  2  and  P  =  3  and  these  differences  are 
now  extended  to  the  upper  surface  of  the  airfoil  where  an  interaction 
between  the  trailing  edge  vortex  and  the  upstream  propagating  shed- 


Under-Resolution  and  Diagnostics  in  Spectral  Simulations 


31 


Figure  21.  Density  contours  of  the  pitching  airfoil  at  non-dimensional  time  t  =  0.75 
corresponding  to  18.55  degrees  angle  of  attack.  Shown  on  the  left  are  contours  at 
spectral  order  P  =  2  and  on  the  right  at  P  =  3. 


Figure  22.  Density  contours  of  the  pitching  airfoil  at  non-dimensional  time  t  =  1.5 
corresponding  to  44.1  degrees  angle  of  attack.  Shown  on  the  left  are  contours  at 
spectral  order  P  =  2  and  on  the  right  at  P  =  3. 
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Figure  23.  Two-dimensional  “z-slice”  of  the  entire  domain  (top  left)  and  detail 
around  the  cylinder  of  the  standard  mesh  (K  =  902  elements  -  top  right);  refined 
(K  =  1,622  elements  -  bottom  left);  and  coarse  mesh  (K  —  412  elements  -  bottom 
right)  used  in  the  spectral  element  simulations.  The  unstructured  grid  shown  is  the 
skeleton  based  on  which  hierarchical  spectral  expansions  are  constructed. 


vortex  takes  place,  as  shown  in  figure  22.  These  flow  pattern  differences 
are  responsible  for  the  aforementioned  differences  in  the  lift  and  drag 
coefficient  at  large  angle  of  attack  as  shown  in  figure  20. 

4.3  Turbulent  Cylinder  Wake 

Numerical  simulation  of  turbulent  wakes  has  been  computationally 
prohibitive  and  only  preliminary  results  have  been  obtained  in  ([8])  us¬ 
ing  DNS.  A  more  systematic  study  of  the  cylinder  turbulent  wake  at 
Re  =  3,900  was  undertaken  by  ([35])  who  used  LES  with  an  upwind 
discretization.  A  second  LES  study  was  performed  by  ([36])  with  central- 
differencing  in  order  to  control  the  numerical  damping  reported  in  the 
first  study,  and  more  recently  a  high-order  LES  study  was  completed 


Under-Resolution  and  Diagnostics  in  Spectral  Simulations 


33 


by  ([37]).  The  results  from  the  three  studies  are  similar  as  far  as  the 
computed  mean  and  rms  velocities  are  concerned,  i.e.  LES  predicts  rel¬ 
atively  accurately,  although  not  uniformly,  the  experimental  results  in 
the  region  downstream  of  x/D  >  3. 

However,  in  the  very-near-wake  all  simulations  converge  to  a  mean 
velocity  profile  in  the  U-shape  unlike  the  experiments  of  ([43])  that  show 
a  V- shape.  In  contrast,  an  independent  LES  study  by  Rodi  and  co¬ 
workers  ([38])  produced  a  V-shape  velocity  profile.  Also,  despite  the 
higher  fluctuations  sustained  in  the  central-differencing  simulations  by  ( 
[39]),  no  clear  inertial  range  was  obtained  in  either  of  the  first  two  LES 
studies  in  contrast  with  the  experiments.  It  is  interesting  to  note  that 
corresponding  simulations  with  the  subfilter  model  turned  off  produced 
an  almost  identical  spectrum  to  the  LES  velocity  spectrum.  A  system¬ 
atic  grid-refinement  study  performed  in  ([35])  also  suggests  that  these 
results  are  resolution-  independent  for  at  least  the  first  ten  diameters  in 
the  near-wake.  The  high-order  LES  of  Kravchenko  &  Moin,  however, 
reproduced  accurately  the  inertial  range  but  predicted  the  same  mean 
velocity  field  (i.e.  U-shape)  as  the  previous  two  simulations. 


Figure  24-  DNS  mean  streamwise  velocity  predictions  at  x/D  =  1.06;  1.54;  2.02 
(from  top  to  bottom,  respectively  ),  (wide  domain  -  solid  line)  and  (narrow  domain  - 
dash  line).  Squares  are  data  of  ([43]). 


For  the  simulations  presented  here  a  continuous  Galerkin  spectral//ip 
element  method  was  employed  in  x-  and  y-directions  while  a  Fourier  ex¬ 
pansion  was  employed  along  the  homogeneous  direction  (cylinder-axis) 
with  appropriate  dealiasing.  Specifically,  triangular  elements  are  used, 
filled  with  Jacobi  polynomial  modes  of  order  P.  We  performed  several 
simulations  corresponding  to  /i-refinement  (i.e.  with  respect  to  number 
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Figure  25.  DNS  rms  streamwise  velocity  predictions  at  x/D  =  1.06;  1.54;  2.02  (from 
top  to  bottom,  respectively  ).  (wide  domain  -  solid  line)  and  (narrow  domain  -  dash 
line).  Squares  are  data  of  ([43]. 


of  elements  K)  and  P-refinement  (with  respect  to  polynomial  order  P ) 
([40]).  In  figure  23  we  show  a  “z-slice”  of  the  computational  domain  in 
the  x  —  y  plane  with  three  different  discretizations.  The  top  plot  shows 
a  grid  with  K  =  902  triangular  prismatic  elements,  which  has  been  the 
standard  grid  we  have  used  for  most  cases.  In  the  bottom  plot  we  also 
show  a  grid  with  finer  resolution  around  the  shear  layers  corresponding 
to  I\  =  1,622  elements,  and  also  a  grid  with  coarser  resolution  cor¬ 
responding  to  K  =  412  elements.  The  polynomial  order  per  element 
varied  from  P  =  4  to  10,  and  the  number  of  Fourier  modes  varied  from 
N  =  2  to  128  (the  corresponding  number  of  physical  points  is  twice  the 
number  of  modes).  The  finest  resolution  simulation  employed  K  =  902 
elements  of  order  P  =  10  and  256  points  {N  =  128  Fourier  modes)  in  the 
spanwise  direction.  The  lowest  resolution  simulation  employed  K  =  412 
elements  with  P  =  6  and  only  N  =  2,  i.e.  a  severe  truncation  of  Fourier 
modes  in  the  spanwise  direction. 

The  domain  extends  from  —15 D  at  the  inflow  to  25 D  at  the  outflow, 
and  from  —9 D  to  9 D  in  the  cross-flow  direction.  Neumann  boundary 
conditions  (i.e.  zero  flux)  were  used  at  the  outflow  and  on  the  sides  of  the 
domain  to  minimize  the  effect  of  normal  boundary  layers  at  the  truncated 
domain.  The  spanwise  length  was  varied  as  Lz/D  =  7r/2,  7t,  1.57T,  27t.  For 
reference,  the  spanwise  length  used  in  all  simulations  of  ([35,  36,  37])  was 
Lz/D  =  ir. 

The  experimental  results  of  ([41])  suggest  a  value  of  correlation  length 
less  than  1.5 D  at  three  diameters  downstream;  this  was  obtained  using 
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the  streamwise  velocity  only.  However,  from  plots  of  the  autocorrelation 
function  for  all  three  velocity  components  and  the  pressure  we  have  seen 
that  at  a  centerline  point  Ruu  drops  to  zero  at  about  1.5 D  but  that, 
in  general,  at  points  off-centerline  Rvv  and  Ruu  do  not  decay  as  fast  ( 
[40]).  Such  results  indicate  that  values  of  Ruu  obtained  in  experiments 
at  centerline  points  may  under-predict  the  spanwise  correlation  length. 
Therefore,  it  may  be  inadequate  to  use  Ruu  as  the  only  criterion  in 
deciding  on  the  domain  size.  Indeed,  we  have  found  that  the  span  length 
is  very  important  in  determining  the  rms  values  in  the  very-near-wake 
and  correspondingly  the  mean  velocity  profiles. 

In  ( [40] )  high  resolution  results  can  be  found  for  many  different  quan¬ 
tities.  Typical  velocity  profiles  for  the  mean  and  the  variance  and  for 
different  spans  are  shown  in  figures  24  and  25.  Here,  we  examine  how 
such  results  are  affected  by  substantially  reducing  the  grid  resolution 
and  without  using  any  subfilter  model.  In  particular,  we  present  here 
results  obtained  on  the  grids  shown  in  figure  23  (bottom  right)  consist¬ 
ing  of  K  =  412  triangular  elements  and  only  P  =  6  and  the  equivalent 
K  =  902  and  P  =  4,  both  cases  corresponding  to  approximately  the 
same  number  of  degrees  of  freedom.  More  comparisons  with  the  subfil¬ 
ter  model  on  the  same  grid  can  be  found  in  ([42]).  We  will  first  use  only 
two  Fourier  modes  in  the  span,  i.e.  the  mean  mode  and  one  perturba¬ 
tion  (IV  =  2  or  4  points).  We  also  choose  a  small  value  for  the  spanwise 
length  Lz/D  =  ir/2. 


Figure  26.  Streamwise  mean  velocity  profile  at  x/D  =  1.06, 1.54,  2.02.  Squares 
denote  experimental  data  of  ([43]),  solid  line  DNS  (K  =  412;  P  =  6),  dash-dot  line 
DNS  (K  =  902;  P  =  4)  and  dash  line  LES  of  Beaudan  &  Moin  ([35]). 

We  compare  first  with  the  experiments  of  ([43])  in  the  very  near- wake 
and  subsequently  with  the  experiments  of  ([41])  farther  downstream. 
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Figure  27.  Mean  velocity  profiles  at  x/D  =  4,7,10.  Circles  denote  experimental 
data  of  ([41]),  solid  line  DNS  (K  =  412;  P  =  6),  dash-dot  line  DNS  (K  =  902;  P  =  4) 
and  dash  line  LES  of  Beaudan  &  Moin  ([35]). 


Figure  28.  Left:  Turbulent  intensity  of  the  streamwise  velocity  (u™s)  at  x/D  = 
4,  7, 10.  Right:  Turbulent  intensity  of  the  cross-flow  velocity  ( »™s )  at  x/D  =  4,  7, 10. 
Circles  denote  experimental  data  of  Ong  &  Wallace,  solid  line  DNS  ( K  =  412;  P  =  6), 
dash-dot  line  DNS  (K  =  902;  P  —  4)  and  dash  line  LES  of  Beaudan  &  Moin  ([35]). 


In  figure  26  we  plot  the  mean  streamwise  velocity  profile  at  locations 
x/D  =  1.06, 1.54  and  2.02.  We  also  include  the  experimental  data  of  ( 
[43])  taken  from  ([35]),  and  the  LES  data  of  ([35]).  We  see  that  the  pre¬ 
dictions  from  both  low-resolution  simulations  without  subfiltering  are 
comparable  to  the  LES  predictions.  In  figures  27,  28  we  plot  the  mean 
streamwise  velocity  and  turbulent  fluctuations,  respectively,  at  locations 
x/D  =  4,  7, 10  and  compare  with  the  experimental  data  of  ([41]).  The 
predictions  for  the  mean  velocities  are  good  but  the  streamwise  turbu¬ 
lence  intensity  shows  some  wiggles,  which  is  an  indication  of  insufficient 
resolution.  However,  the  low-resolution  spectral  simulations  obtain  an 
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overall  better  agreement  with  the  experimental  data  than  the  dissipative 
LES  predictions  reported  in  ([35]). 


Figure  29.  Streamwise  mean  velocity  profile  at  x/D  =  1.06  for  different  Fourier 
modes  employed  along  the  span.  Squares  denote  experimental  data  of  ([43]),  dash- 
dot  line  2D  simulation,  dash  line  N  =  2,  dot-solid  line  TV  =  8  and  solid  line  TV  =  32 
(coincides  with  N  =  8). 


The  results  presented  so  far  were  obtained  with  only  N  =  2  Fourier 
modes  employed  along  the  cylinder  span.  Of  interest  is  to  examine  the 
influence  of  the  number  of  Fourier  modes  N  on  the  mean  velocity  profiles 
presented  above  while  retaining  the  same  resolution  in  the  x  —  y  planes. 
We  performed  additional  simulations  with  N  =  8  and  32  and  also  a  two- 
dimensional  simulation.  As  we  see  in  figure  29  there  is  essentially  no 
difference  in  the  predicted  mean  streamwise  velocity  profile  from  N  =  2 
to  N  =  2>2  but  the  two-dimensional  prediction  deviates  substantially. 
The  cases  with  N  =  8  and  N  =  32  correspond  to  almost  identical 
predictions  suggesting  convergence  in  the  z-direction. 

The  results  presented  here  indicate  that  the  first  Fourier  mode  carries 
most  of  the  spanwise  energy  for  the  chosen  span  Lz/D  =  ir/2,  as  it  is 
evident  by  comparing  with  the  two-dimensional  results  in  figure  29.  This 
has  been  independently  verified  by  computing  the  averaged  plane-modal 
energy  Exy(m)  =  /  [u^  +  v2^  +  w^]  dxdy  and  observe  its  decay  with 
respect  to  the  mode  number. 

Given  the  surprisingly  good  results  with  this  low-resolution  at  Re  = 
3, 900,  we  performed  another  set  of  simulations  with  the  same  low- 
resolution  and  with  only  N  =  2  Fourier  modes  at  Re  =  5,  000  for  which 
we  had  available  experimental  data  from  the  work  of  ([44]).  In  figure  30 
we  plot  the  mean  velocity  profile  and  the  streamwise  turbulent  intensity 
u'2  at  station  x/D  =  10.  Again,  we  see  that  despite  some  wiggles  in  the 
numerical  results  the  agreement  with  the  experimental  results  is  good. 
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Figure  30.  Streamwise  mean  velocity  profile  (bottom)  and  turbulent  fluctuation 
(top)  at  x/D  =  10  and  Re  =  5,000.  The  experimental  data  (circles)  are  from  Zhou 
&  Antonia  (1993). 


5.  Discussion 

Spectral  methods  have  been  used  with  great  success  in  simulating 
turbulent  flows  in  periodic  cubes  and  channel  domains.  The  algorith¬ 
mic  developments  of  the  last  two  decades  have  led  to  a  new  simulation 
capability  for  complex-geometry  turbulent  flows  as  well.  Such  capabil¬ 
ity,  in  conjunction  with  terascale  computing  at  the  PC  cluster  level, 
will  undoubtly  lead  to  significant  advances  in  simulating  turbulence  in 
more  realistic  configurations  and  in  realistic  operating  conditions.  In 
this  chapter,  we  have  summarized  some  of  these  developments  and  have 
presented  results  for  several  transitional  and  turbulent  flows  in  complex- 
geometry  domains. 

Great  care  has  to  be  exercised  however  in  interpreting  results  from 
DNS  or  LES  in  simple-geometry  flows  at  high  Reynolds  number  or  in 
complex-geometry  flows  with  new  physics.  Many  of  these  simulations 
may  be  under-resolved  at  some  level;  for  example,  although  the  velocity 
mean  and  variance  may  be  correctly  predicted,  the  high-order  statistics 
or  the  dissipation  spectrum  may  be  erroneous.  It  is  important,  therefore, 
to  have  a  diverse  set  of  diagnostic  tools  to  characterize  numerical  uncer- 
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tainty  in  these  situations.  In  particular,  understanding  how  numerical 
methods  behave  in  DNS  and  LES  for  many  prototype  flows  provides  an 
insight  into  uncert amities  and  their  origin  in  large-scale  simulations.  Val¬ 
idation  and  verification  ([1])  of  a  turbulent  simulation  is  a  very  difficult 
task  that,  unfortunately,  cannot  be  handled  solely  by  error  estimators 
which  are  based  primarily  on  extensions  of  linear  concepts  ([45,  46,  47]). 

Under-resolved  simulations  are  not  useless  if  the  numerical  uncertainty 
is  properly  characterized  ([48,  49,  50]),  i.e.  quantified  with  a  properly 
constructed  error  bar  ([51]).  In  many  cases  such  under-resolved  simula¬ 
tions  may  contain  the  answer  that  we  seek,  e.g.  an  averaged  lift  or  drag 
coefficient  and  at  the  accuracy  level  that  we  expect.  The  type  of  flow 
that  we  simulate  is  important  in  that  respect.  We  have  presented  here, 
for  example,  low-resolution  simulations  of  the  cylinder  turbulent  wake 
which  lead  to  results  that  match  the  experiments  at  Re  =  5,  000.  Flows 
with  inherent  low-dimensionality,  such  as  the  cylinder  wake  for  which  the 
vortex  shedding  process  dominates  the  dynamics,  may  then  be  easier  to 
simulate.  Moreover,  quantification  of  numerical  uncertainty  in  hierar¬ 
chical  methods,  such  as  the  spectral//ip  element  method,  is  also  easier 
to  achieve  and  offer  the  possibility  of  obtaining  multiple  solutions  with 
relatively  simple  P-refinement  without  the  need  for  re-meshing  which, 
typically,  is  a  large  overhead  component  in  a  turbulence  simulation. 
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